{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem Set 3 - Finding planets outside the solar system" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Total points:** 30\n", "\n", "**Due:** Friday 18th September 2020, 7pm CEST\n", "\n", "**Format:** Jupyter Notebook or python program" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Background" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this problem set, we are going to investigate a dataset that contains the *radial velocity* of a star at various points in time. The radial velocity is the velocity with which the star is moving towards or away from us (positive velocities indicate it is moving away from us).\n", "\n", "If a star is not close to any other objects, there is no reason why the radial velocity should change over time. However, in the case where a second object is orbiting a star (such as a planet or another star), the star and the object will both orbit the center of mass of the system. Therefore, the star will show periodic variations in its velocity over time. These changes in velocity then cause a shift in spectral lines via the Doppler effect, which we can measure with telescopes equipped with state-of-the-art spectrographs. The smaller the second object, the less the star will be affected. For example, the Earth causes the Sun to change its velocity with an amplitude of 0.1 m/s over 1 year.\n", "\n", "See the following video to see an example of a large planet orbiting a star and the effect on the observed spectral lines of the star:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.display import YouTubeVideo\n", "YouTubeVideo(\"-BuwWtMygxU\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Measuring the radial velocities of stars to very high precision can therefore be used to find planets. Note that we cannot simply measure the velocity of planets around other stars, since they are much fainter than the stars themselves.\n", "\n", "In this problem set, we want to find out whether a particular star has a companion object, and if so, we want to estimate the probability that the companion object is a planet. The data file required for this exercise is ``data/UID_0113357_RVC_001.tbl`` from the ZIP-file from the course homepage or ``data/PS03-data/UID_0113357_RVC_001.tbl.txt`` on the KIP Jupyter Server. Pick the file of your choice and then carry out the analysis described in Parts 1 and 2.\n", "\n", "**The first column of the file contains the time in days, the second column contains the radial velocity in m/s, and the third column contains the uncertainty in the radial velocity in m/s.** The fourth column can be ignored." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 1 (15 points)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start off, we want to see whether the star does indeed show periodic variations, and if so, we want to measure the period and amplitude.\n", "\n", "**Read in the data and make a plot of the radial velocities over time** (make sure you include axis labels, units, and error bars!).\n", "\n", "As you will see, the observations were not taken with equal spacing in time, and it does not look like we can easily see any periodic variations by eye, so we will instead use an automated method. The best way to find the period of a time series like this is to use a periodogram algorithm. One algorithm that works well in cases like this is the following:\n", "\n", "* pick a ``period``\n", "\n", "* compute the phase of the radial velocity curve:\n", "\n", " phase = (time % period) / period\n", "\n", "* sort the phase and velocity points by the phase\n", "\n", "* find the \"path-length\" of a line that would join all the points, from left to right. If we write the sorted phase as $\\phi$, the velocity as $v$, and the number of points as $N$, then this is:\n", "\n", " $$\\sum_i^{N-1}{\\sqrt{(\\phi_{i+1} - \\phi_{i})^2 + (v_{i+1} - v_{i})^2}}$$\n", "\n", "* repeat for different periods, and minimize the string length\n", "\n", "**Implement this algorithm** and find the *path length* for 10,000 periods logarithmically spaced between 1 and 100 days. Note that the calculation of the path length should not use any loops, otherwise it will be very slow (but you can still loop over the periods).\n", "\n", "**Hint 1**: you can find ``phi[i+1] - phi[i]`` for all elements of an array with ``np.diff(phi)`` or remember how we computed similar quantities previously.\n", "\n", "**Hint 2**: if you want to sort an array ``y`` according to the values of ``x``, you can use [np.argsort](http://docs.scipy.org/doc/numpy/reference/generated/numpy.argsort.html):\n", "\n", " y_new = y[np.argsort(x)]\n", "\n", "**Make a log-log plot of the string length versus the period** and describe the plot. You should see several clear minima. **Find the period** of the first clear minimum (optional question: why the first? what do the other minima correspond to?).\n", "\n", "**Make a plot** of radial velocity (with error bars) versus phase for the period you found above, and if you have picked the period correctly you should see something that looks similar to one period of a sinusoidal function. This means that there are indeed periodic (and in fact sinusoidal) variations in the radial velocity curve, indicating that a second object is likely orbiting the star!\n", "\n", "**Fit the radial velocity versus phase** using a function of the form:\n", "\n", "$$f(x) = a \\sin{(2\\pi x + b)}$$\n", "\n", "and **make a plot** of the radial velocity (with error bars) and with the best-fit overplotted.\n", "\n", "The parameter $a$ will give the amplitude of the radial velocity curve (if it gives a negative value, be sure to make it positive for the next section). **Find the uncertainty on $a$** using the covariance matrix returned by the fitting function.\n", "\n", "You should now have the period and the amplitude (with uncertainty) of the radial velocity curve!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2 (15 points)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this section, we will now try and estimate the mass of the invisible companion that is orbiting the star. At the end of Part 1, we found that the radial velocity curve could be well fitted by a sinusoidal curve, which suggests that the invisible companion is on a circular orbit (at least to a good approximation). If you were unable to solve Part 1, assume a period of $P=4.5\\,\\mathrm{d}$ and an RV amplitude of $K=60\\pm 1\\,\\mathrm{m}\\,\\mathrm{s}^{-1}$ in the following.\n", "\n", "For two bodies (in our case a star and an unknown object) orbiting a common center of mass in circular orbits, Kepler's 3rd law of motion states that:\n", "\n", "$$\\left(\\frac{2\\pi}{P}\\right)^2 = \\frac{G(M_\\star + m_\\mathrm{P})}{a^3} = \\frac{G M_\\star (1 + q)}{a^3}$$\n", "\n", "where $P$ is the orbital period, $M_\\star$ is the mass of the star, $q$ is the ratio of the object mass to the mass of the star ($q=m_\\mathrm{P}/M_\\star$), and $a$ is the distance between the object and the star (the orbital separation).\n", "\n", "In addition, for two objects orbiting a common center of mass, the radial-velocity amplitudes $K$ are related to the mass of the object such that the ratio is:\n", "\n", "$$\\frac{K_\\star}{K_{\\rm object}} = q$$\n", "\n", "Finally, for a circular orbit (i.e. a constant velocity along the orbit):\n", "\n", "$$K_{\\rm object} P = 2\\pi a$$\n", "\n", "Substituting this back into Kepler's law and simplifying gives:\n", "\n", "$$(1 + q)~q^3 = \\frac{P K_\\star^3}{2\\pi G M_\\star}$$\n", "\n", "Finally, there is one important thing to consider, which is that the velocity amplitude we measure for the star is not $K_\\star$ but\n", "\n", "$$K_{\\star, \\mathrm{obs}} = K_\\star~\\sin{i}$$\n", "\n", "where $i$ is the viewing/inclination angle of the system ($i=90^\\circ$ means that we are observing the system *edge-on*, and $i=0^\\circ$ corresponds to a *face-on* orientation where we would not see the star move towards or away from us). $K_{\\star, \\mathrm{obs}}$ is the amplitude of the velocity - the ``a`` parameter - that you measured in Part 1.\n", "\n", "Therefore, assuming $q << 1$, the final equation is:\n", "\n", "$$q \\approx \\left(\\frac{P}{2\\pi G M_\\star}\\right)^{1/3} \\frac{K_{\\star, \\mathrm{obs}}}{\\sin{i}}$$\n", "\n", "**Compute the value of $q$** assuming the values of $P$ and $K_{\\star, \\mathrm{obs}}$ that you found in Part 1, and assuming $\\sin{i}=1$ and that $M_\\star$ is the mass of the Sun ($M_\\star=1.989\\times 10^{30}$kg). Using $q$, derive the mass of the object in units of the mass of Jupiter ($M_\\mathrm{Jupiter}=1.898\\times 10^{27}$kg). Is the object likely to be a planet? (*Hint:* another very useful package is the [SciPy constants package](https://docs.scipy.org/doc/scipy/reference/constants.html) that contains numerical values of many natural and derived constants that are useful for scientists in general. For those ineterested in Astrophysics, you may also want to check out the [AstroPy constants](http://docs.astropy.org/en/stable/constants/) package.)\n", "\n", "In reality, your measurement of $K_{\\star, \\mathrm{obs}}$ has uncertainties, the mass of the star is uncertain, and the viewing angle is also uncertain. **Carry out a Monte-Carlo error propagation** simulation to find the likely distribution of masses for the object assuming that:\n", "\n", "* the error on the radial velocity is the one derived in Part 1 and is a *normal* error (i.e. the distribution follows a Gaussian)\n", "\n", "* the mass of the star is sampled from a uniform distribution between 0.6 and 1.4 times the mass of the Sun.\n", "\n", "* the viewing angle can be anywhere between 0 and 90 degrees. However, one cannot simply sample $i$ randomly between 0 and 90 because some viewing angles are more likely than others (do you see why?) - to do it properly you need to sample $\\cos{i}$ uniformly between 0 and 1, then derive $i$ or $\\sin{i}$ from this (i.e. the probability density function of the inclination angle is $\\sin i$).\n", "\n", "From the Monte-Carlo simulation, **plot a histogram** of the probability that the object has a certain mass, and show only the range from 0 to 13 times the mass of Jupiter.\n", "\n", "**What is the (numerical) probability that the object is less massive than 13 times the mass of Jupiter?** (this is usually considered to be the upper limit of the planet mass). What degree of confidence do we have that the object is a planet, using the 1/2/3/4/5-sigma confidence terminology? (see [here](http://en.wikipedia.org/wiki/68%E2%80%9395%E2%80%9399.7_rule) for more information). *Hint:* You will have to convert from probability to $\\sigma$ using a Gaussian distribution with $\\mu=0$ and $\\sigma=1$. How can you best do this? We already discussed the required tools in our classes but, alternatively, maybe this [SciPy package](https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.stats.norm.html) comes in handy... \n", "\n", "Based on this, what can you conclude about the object?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.4" } }, "nbformat": 4, "nbformat_minor": 4 }